home *** CD-ROM | disk | FTP | other *** search
/ Power Bytes: Money & Finance / PowerBytes Money and Finance CD-ROM 01 / PowerBytes Money and Finance CD-ROM 01.iso / HyperCard Files / MacMathPascal / card_8956.txt < prev    next >
Encoding:
Text File  |  1988-09-05  |  7.1 KB  |  315 lines
  1. -- card: 8956 from stack: in
  2. -- bmap block id: 0
  3. -- flags: 0000
  4. -- background id: 2607
  5. -- name: EIGEN.PAS
  6.  
  7.  
  8. -- part contents for background part 22
  9. ----- text -----
  10. EIGEN.PAS
  11.  
  12. -- part contents for background part 23
  13. ----- text -----
  14. EIGEN.PAS
  15.  
  16. -- part contents for background part 26
  17. ----- text -----
  18. 14
  19.  
  20. -- part contents for background part 17
  21. ----- text -----
  22. EIGEN  [procedure]
  23.  
  24.  
  25.  
  26. PURPOSE
  27.  
  28. EIGEN finds eigenvalues and eignevectors for the eigenvalues from a matrix "AI".   An iteration procedure is used to find the smallest or largest eigenvalue and then the associated eigenvector.   Then "S" other eigenvalues are found sweeping.  (S is an input parameter.)  There is a boolean input parameter "Small" which determines if the smallest or largest eigenvalue is found first.  If Small is TRUE then the smallest eigenvalue is found followed by the rest of the "S" eigenvalues.  If small is TRUE the the matrix AI is inverted as a by product of the procedure to find the eigenvalues.
  29.  
  30.  
  31. TYPE REQUIREMENTS
  32. CONST
  33. MAXSIZE = THE MAXIMUM SIZE OF THE INPUT ARRAY
  34. MAXPLUS1 = MAXSIZE +1
  35.  
  36. Type
  37.   FLOAT = extended or real;
  38.   FIXXED = INTEGER;
  39.   SQUARE = array[1..MAXSIZE, 1..MAXSIZE] of FLOAT;
  40.   MYRECT = array[1..MAXPLUS1, 1..MAXSIZE] of FLOAT;
  41.   RAY = array[1..MAXSIZE] of FLOAT;
  42.  
  43.  
  44. CALLING PROCEDURE
  45.  
  46. var
  47. AI:SQUARE;{matrix of interest}
  48. Z:Square;{output matrix of eigenvectors};
  49. Lambda:ray;{Array of eigenvalues}
  50. N:fixxed;{size of the matrix}
  51. S:FIXXED;{How many eigenvalues to find}
  52. Epsi:float;{Accuracy requirement.  determines when to stop iteration in finding the first eigenvalues}
  53. SMALL:boolean; {if small is true the smallest eigenvalue will be found first and AI will be inverted.  If small is false the the largest eigenvalue will be found first}
  54. ITERATE:FIXXED;{Number of times MULTI_DIFEQ will iterate before returning}
  55.  
  56. Define the procedure Eigen.  The calling sequence to Eigen is 
  57.  
  58.  EIGEN(AI, LAMBDA, Z, N, S, EPSI, SMALL,ITERATE);
  59.  
  60. Define AI,N,S,EPSI and SMALL.   Then call the procedure.
  61.  
  62.  EIGEN(AI, LAMBDA, Z, N, S, EPSI, SMALL,ITERATE);
  63.  
  64. When the procerdure finishes, the matrix "Z" will contain the eigenvectors and "Lambda" will contain the eigenvalues.  The N components of the eigenvector of corresponding to the first eigenvalue are in Z[1,1],Z[1,2]..Z[1,N].
  65.  
  66.  
  67. EXAMPLE
  68.  
  69. The example defines a 3 x 3 matrix and finds 3 eigenvalues and eigenvectors starting with the largest.  The eigenvalues and vectors are printed.
  70.  
  71.  
  72. REFERENCES
  73. James, M. L., Smith, G. M., Wolford, J. C.:  Applied Numerical Methods for 
  74. Digital Computation With Fortran and CSMP, Harper and Row, New York, 
  75. 1977.
  76.  
  77.  
  78. -- part contents for background part 6
  79. ----- text -----
  80. PROGRAM DOEIGEN;
  81. CONST
  82. MAXSIZE = 10;
  83. MAXPLUS1 = 11;
  84. TYPE
  85. FLOAT = extended;
  86. FIXXED = INTEGER;
  87. SQUARE = ARRAY[1..MAXSIZE, 1..MAXSIZE] OF FLOAT;
  88. MYRECT = ARRAY[1..MAXPLUS1, 1..MAXSIZE] OF FLOAT;
  89. data = ARRAY[1..MAXSIZE] OF FLOAT;
  90. VAR
  91. inmatrix, eigenvector : SQUARE;
  92. N, S, SP1, I, L, ITMAX : FIXXED;
  93. SMALL : BOOLEAN;
  94. EPSI : FLOAT;
  95. LAMBDA : data;
  96. PROCEDURE see;
  97. VAR
  98. R : Rect;
  99. BEGIN
  100. HideAll;
  101. SetRect(R, 0, 35, 550, 330);
  102. SettextRect(R);
  103. Showtext;
  104. END;
  105. PROCEDURE eigen (VAR inmatrix : SQUARE;
  106. VAR LAMBDA : data;
  107. VAR eigenvector : SQUARE;
  108. N, SIN : FIXXED;
  109. EPSI : FLOAT;
  110. OVER : BOOLEAN;
  111. ITMAX : FIXXED);
  112. VAR
  113. BACK : ARRAY[1..maxsize, 1..2] OF FIXXED;
  114. K, JJ, KP1, I, J, L, KROW, IROW : FIXXED;
  115. rot, TEMP : FLOAT;
  116. S : INTEGER;
  117. X, D, C, G, ADD, Y, XX : data;
  118. REC1, REC2, REC3 : MYRECT;
  119. SUM : FLOAT;
  120. ITERATE, IP1, II, KK : FIXXED;
  121. done10, done20, done30, done40 : BOOLEAN;
  122. BEGIN
  123. done10 := FALSE;
  124. done20 := FALSE;
  125. done30 := FALSE;
  126. done40 := FALSE;
  127. S := SIN - 1;
  128. IF (OVER) THEN
  129. BEGIN
  130. FOR K := 1 TO n DO
  131. BEGIN
  132. JJ := K;
  133. IF (K <> n) THEN
  134. BEGIN
  135. KP1 := K + 1;
  136. rot := ABS(inmatrix[K, K]);
  137. FOR I := KP1 TO n DO
  138. BEGIN
  139. TEMP := ABS(inmatrix[I, K]);
  140. IF (rot < TEMP) THEN
  141. BEGIN
  142. rot := TEMP;
  143. JJ := I;
  144. END;
  145. END;
  146. END;
  147. BACK[K, 1] := K;
  148. BACK[K, 2] := JJ;
  149. IF (JJ <> K) THEN
  150. FOR J := 1 TO n DO
  151. BEGIN
  152. TEMP := inmatrix[JJ, J];
  153. inmatrix[JJ, J] := inmatrix[K, J];
  154. inmatrix[K, J] := TEMP;
  155. END;
  156. FOR J := 1 TO n DO
  157. IF (J <> K) THEN
  158. inmatrix[K, J] := inmatrix[K, J] / inmatrix[K, K];
  159. inmatrix[K, K] := 1.0 / inmatrix[K, K];
  160. FOR I := 1 TO n DO
  161. IF (I <> K) THEN
  162. FOR J := 1 TO n DO
  163. IF (J <> K) THEN
  164. inmatrix[I, J] := inmatrix[I, J] - inmatrix[K, J] * inmatrix[I, K];
  165. FOR I := 1 TO n DO
  166. IF (I <> K) THEN
  167. inmatrix[I, K] := -inmatrix[I, K] * inmatrix[K, K];
  168. END;
  169. FOR L := 1 TO n DO
  170. BEGIN
  171. K := n - L + 1;
  172. KROW := BACK[K, 1];
  173. IROW := BACK[K, 2];
  174. IF (KROW <> IROW) THEN
  175. FOR I := 1 TO n DO
  176. BEGIN
  177. TEMP := inmatrix[I, KROW];
  178. inmatrix[I, KROW] := inmatrix[I, IROW];
  179. inmatrix[I, IROW] := TEMP;
  180. END;
  181. END;
  182. END;
  183. FOR I := 1 TO N DO
  184. BEGIN
  185. X[I] := 1.0 / SQR(1.0 * N) + 1.0 / (1.0 * I);
  186. END;
  187. ITERATE := 0;
  188.  
  189. REPEAT
  190. FOR I := 1 TO N DO
  191. BEGIN
  192. D[I] := 0.0;
  193. FOR J := 1 TO N DO
  194. D[I] := D[I] + inmatrix[I, J] * X[J];
  195. END;
  196. ITERATE := ITERATE + 1;
  197. FOR I := 1 TO N DO
  198. eigenvector[1, I] := D[I] / D[1];
  199. I := 1;
  200. REPEAT
  201.  
  202. IF (ABS(X[I] - eigenvector[1, I]) - EPSI * eigenvector[1, I] >= 0) THEN
  203. done10 := TRUE;
  204. I := I + 1;
  205. UNTIL (done10 OR (I > N));
  206. IF (NOT (done10)) THEN
  207. done20 := TRUE;
  208.  
  209. IF (NOT (done20)) THEN
  210. BEGIN
  211. FOR I := 1 TO N DO
  212. X[I] := eigenvector[1, I]
  213. END;
  214. IF (ITERATE > ITMAX) THEN
  215. done20 := TRUE;
  216. UNTIL (done20);
  217.  
  218. LAMBDA[1] := 1.0 / D[1];
  219. FOR I := 1 TO S DO   {100 START}
  220. BEGIN
  221. IP1 := I + 1;
  222. FOR L := 1 TO N DO
  223. REC1[IP1, L] := 1.0;
  224. G[I] := 0.0;
  225. FOR L := 1 TO N DO
  226. G[I] := G[I] + eigenvector[I, L] * eigenvector[I, L];
  227. ITERATE := 0;
  228. REPEAT
  229. done30 := FALSE;
  230.  
  231. FOR II := 1 TO I DO
  232. BEGIN
  233. SUM := 0.0;
  234. FOR L := 1 TO N DO
  235. SUM := SUM + eigenvector[II, L] * REC1[IP1, L];
  236. C[II] := SUM / G[II];
  237. END;
  238. ITERATE := ITERATE + 1;
  239. FOR L := 1 TO N DO
  240. ADD[L] := 0.0;
  241. FOR K := 1 TO I DO
  242. BEGIN
  243. FOR L := 1 TO N DO
  244. BEGIN
  245. Y[L] := C[K] * eigenvector[K, L];
  246. ADD[L] := ADD[L] + Y[L];
  247. END;
  248. END;
  249. FOR L := 1 TO N DO
  250. REC2[IP1, L] := REC1[IP1, L] - ADD[L];
  251. FOR II := 1 TO N DO
  252. BEGIN
  253. XX[II] := 0;
  254. FOR KK := 1 TO N DO
  255. XX[II] := XX[II] + inmatrix[II, KK] * REC2[IP1, KK];
  256. END;
  257. FOR L := 1 TO N DO
  258. REC3[IP1, L] := XX[L] / XX[1];
  259. L := 1;
  260. done40 := FALSE;
  261. REPEAT
  262. IF ((ABS(REC1[IP1, L] - REC3[IP1, L]) - EPSI * REC3[IP1, L]) > 0.0) THEN
  263. done40 := TRUE;
  264. L := L + 1;
  265. UNTIL ((L > N) OR (done40));
  266. IF (done40) THEN
  267. BEGIN  {24}
  268. FOR L := 1 TO N DO
  269. REC1[IP1, L] := REC3[IP1, L];
  270. END;
  271. IF (ITERATE > ITMAX) THEN
  272. done30 := TRUE;
  273. UNTIL (done30);
  274. IF (over) THEN
  275. LAMBDA[IP1] := 1.0 / XX[1]
  276. ELSE
  277. LAMBDA[IP1] := XX[1];
  278.  
  279. FOR L := 1 TO N DO
  280. eigenvector[IP1, L] := REC3[IP1, L];
  281. END;
  282. IF (NOT over) THEN
  283. lambda[1] := 1.0 / lambda[1];
  284. END;
  285. BEGIN
  286. see;
  287. inmatrix[1, 1] := 20.0;
  288. inmatrix[1, 2] := -10.0;
  289. inmatrix[1, 3] := 0.0;
  290. inmatrix[2, 1] := -10.0;
  291. inmatrix[2, 2] := 20.0;
  292. inmatrix[2, 3] := -10.0;
  293. inmatrix[3, 1] := 0.0;
  294. inmatrix[3, 2] := -10.0;
  295. inmatrix[3, 3] := 10.0;
  296. N := 3;
  297. S := 3;
  298. EPSI := 1E-3;
  299. ITMAX := 50;
  300. writeln(' input matrix');
  301. writeln(inmatrix[1, 1] : 10 : 5, inmatrix[1, 2] : 10 : 5, inmatrix[1, 3] : 10 : 5);
  302. writeln(inmatrix[2, 1] : 10 : 5, inmatrix[2, 2] : 10 : 5, inmatrix[2, 3] : 10 : 5);
  303. writeln(inmatrix[3, 1] : 10 : 5, inmatrix[3, 2] : 10 : 5, inmatrix[3, 3] : 10 : 5);
  304. writeln;
  305. SMALL := false;
  306. EIGEN(inmatrix, LAMBDA, eigenvector, N, S, EPSI, SMALL, ITMAX);
  307. FOR I := 1 TO S DO
  308. BEGIN
  309. WRITELN;
  310. WRITELN('  LAMBDA(', I : 2, ')=', LAMBDA[I] : 20 : 10);
  311. WRITELN(' THE  EIGENVECTOR IS ');
  312. FOR L := 1 TO N DO
  313. WRITELN(eigenvector[I, L] : 20 : 10);
  314. END;
  315. END.